library(clusterProfiler)
simpl_enrich <- simplify
source('00_util_scripts/mod_bplot.R')

proj.nm <- 'mission/FPP/xiangya_sle_scRNA/'

# examine mva pathway ----------
kegg_mva <-
  c('Acat1','Acat2','Hmgcs1','Hmgcs2','Hmgcr','Mvk','Pmvk','Mvd','Idi1','Idi2','Fdps') |>
  str_to_upper()

# examine interested cytokines ------
key_cytokine <- c('Il6','Ccl20','Tslp','Flt3lg','Csf2','Tnf') |>
  str_to_upper()

kera_slevhc <- read_csv('mission/FPP/xiangya_sle_scRNA/hs_kera_sle-hc_deg.csv')

sle_v3hvl <- read_csv('mission/FPP/xiangya_sle_scRNA/hs_sle_v3hvl.csv')

hc.v3hvl <- read_csv('mission/FPP/xiangya_sle_scRNA/hs_hc_v3hvl.deg.csv')

v3h_slevhc <- read_csv('mission/FPP/xiangya_sle_scRNA/hs_v3h_sle-hc_deg.csv')

v3l_slevhc <- read_csv('mission/FPP/xiangya_sle_scRNA/hs_v3l_sle-hc_deg.csv')

fine.type.slevhc.deg <-
  read_csv('mission/FPP/xiangya_sle_scRNA/results/mix4.fine.slevhc.deg.csv')

# KC: SLE vs HC ------------
### FIG: volcano ============
kera_slevhc |>
  mutate(p_val_adj = ifelse(gene == 'CCL20', 1e-295, p_val_adj)) |>
  filter(p_val_adj != 0) |>
  plot_pub_volc(group1 = 'SLE KC', group2 = 'HC KC',
                highlights = c('CCL20','TNF','FLT3LG','TRPV3'), force = TRUE)

publish_pdf('mission/FPP/figures/SLE-HC_kera_volcano.pdf', width = 65)

kera_slevhc_cyto <- kera_slevhc |>
  filter(gene %in% key_cytokine)

### FIG: up ORA go =======
kera_sle_upgo <- fine.type.slevhc.deg |>
  filter(avg_log2FC > 1 & p_val_adj < .05, cluster == 'KC') |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Hs.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           qvalueCutoff = .05,
           readable = T) |>
  simpl_enrich()

kera_sle_upgo@result |>
  write_csv('mission/FPP/xiangya_sle_scRNA/kc.slevhc.upgo.csv')

kera_sle_upgo@result |>
  publish_enrichment()

publish_pdf('mission/FPP/figures/mix4.epi.kc.SLE.vs.HC.upgo.ora.pdf',
            width = 55)

kera_slevhc |>
  right_join(chole_p_tib, join_by(gene == chole_p_str))

### GSEA -------
kc.slevhc.gsego <- kera_slevhc |>
  filter(p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(OrgDb = 'org.Hs.eg.db',
        keyType = 'SYMBOL',
        ont = 'BP')

kc.slevhc.gsego@result |>
  filter(str_detect(Description, 'oxygen|stress')) |>
  as_tibble()

kc.slevhc.gsego |>
  gseaplot2('GO:0072593',
            title = 'GO pathway: reactive oxygen species metabolic process')

# HC: v3h vs v3l -----------
### FIG: volcano ==========
hc.v3hvl |>
  filter(p_val_adj != 0) |>
  plot_bill_volc('HC V3-hi KC', 'HC V3-lo KC',
                 highlights = c(key_cytokine, kegg_mva)) +
  ggtitle('HC TRPV3-hi KC vs HC TRPV3-lo KC')

hc.v3hvl |>
  filter(p_val_adj != 0) |>
  plot_pub_volc('HC V3-hi KC', 'HC V3-lo KC',
                 highlights = c(key_cytokine, kegg_mva))

publish_pdf('mission/FPP/figures/HC_v3h-v3l_volcano.pdf')

### FIG: up ORA go ==========
hc.v3hvl_upgo <- hc.v3hvl |>
  filter(avg_log2FC > 1 & p_val_adj < .05) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Hs.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T,
           qvalueCutoff = .05) |>
  simpl_enrich()

hc.v3hvl_upgo@result |>
  plot_enrichment() +
  ggtitle('Upregulated GO BP pathway in HC V3-hi KC vs HC V3-lo KC')

hc.v3hvl_upgo@result |>
  head(10) |>
  publish_enrichment() +
  ggtitle('Upregulated GO BP pathway in HC V3-hi KC vs HC V3-lo KC')

publish_pdf('mission/FPP/figures/HC_v3h-v3l_upgo.pdf', width = 65)

hc.v3hvl_upgo <- hc.v3hvl |>
  filter(avg_log2FC > 1 & p_val_adj < .05) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Hs.eg.db',
           keyType = 'SYMBOL',
           ont = 'ALL',
           readable = T,
           qvalueCutoff = .05) |>
  simpl_enrich()

hc.v3hvl_upgo@result |>
  write_csv('mission/FPP/xiangya_sle_scRNA/hc.v3hvl.upgo.csv')

read_tsv('mission/FPP/xiangya_sle_scRNA/hc.v3hvl.go.digest.tsv') |>
  plot_enrichment(n = 12) +
  labs(title = 'Upregulated GO pathways in human HC v3h vs v3l KC')

### up ORA KEGG =========
hc.v3hvl_upkegg <- hc.v3hvl |>
  filter(avg_log2FC > 0 & p_val_adj < .05) |>
  pull(gene) |>
  bitr(fromType = 'SYMBOL', toType = 'ENTREZID',
       OrgDb = 'org.Hs.eg.db') |>
  pull(ENTREZID) |>
  enrichKEGG(qvalueCutoff = .1)

hc.v3hvl_upkegg@result |>
  filter(qvalue < .1) |>
  write_csv('mission/FPP/xiangya_sle_scRNA/hc.v3hvl.upkegg.csv')

hc.v3hvl_upkegg@result |>
  filter(qvalue < .1) |>
  plot_enrichment() +
  labs(title = 'Upregulated KEGG pathways in human HC v3h vs v3l KC')

### GO GSEA ===========
hc.v3hvl_gse.kegg <- hc.v3hvl |>
  pull(gene) |>
  bitr(fromType = 'SYMBOL', toType = 'ENTREZID',
       OrgDb = 'org.Hs.eg.db') |>
  left_join(hc.v3hvl, join_by(SYMBOL == gene)) |>
  filter(p_val_adj < .05) |>
  pull(avg_log2FC, name = ENTREZID) |>
  sort(decreasing = T) |>
  gseKEGG()

hc.v3hvl_gse.kegg@result |>
  DT::datatable()

# SLE: v3h vs v3l -----------
### FIG: volcano =========
sle_v3hvl |>
  plot_bill_volc('SLE V3-hi KC', 'SLE V3-lo KC',
                highlights = c('HMGCR','TRPV3'), force = T)

sle_v3hvl |>
  mutate(p_val_adj = ifelse(gene == 'TRPV3', 1e-300, p_val_adj)) |>
  filter(abs(avg_log2FC) > .1, p_val_adj != 0) |>
  plot_pub_volc('SLE V3-hi KC', 'SLE V3-lo KC',
                highlights = c('CCL20', 'TRPV3', 'HMGCR'), force = T)

g1 <- last_plot()

sle_v3hvl.up5 <- g1$data |>
  filter(p_val_adj < 1e-80, avg_log2FC > 6) |>
  mutate(gene = ifelse(str_detect(gene, 'BX'), 'ZNG1F', gene)) |>
  slice_max(avg_log2FC, n = 3, with_ties = F)

sle_v3hvl.sig5 <- g1$data |>
  filter(p_val_adj < 1e-100, avg_log2FC < -2.5) |>
  slice_min(p_val_adj, n = 4, with_ties = F) |>
  bind_rows(sle_v3hvl.up5)

g1 +
  geom_point(data = sle_v3hvl.sig5, size = 1, color = 'darkgreen') +
  geom_text_repel(data = sle_v3hvl.sig5, aes(label = gene), color = 'black',
                  size = 2, box.padding = .1)

publish_source_plot('SLE_v3h-v3l_volc_v2')

### FIG: up ORA go ==========
sle_v3hvl_upgo <- sle_v3hvl |>
  filter(avg_log2FC > 1 & p_val_adj < .05) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Hs.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T,
           qvalueCutoff = .05) |>
  simpl_enrich()

sle_v3hvl_upgo@result |>
  publish_enrichment()

publish_pdf('mission/FPP/figures/mix4.epi.SLE_v3h-v3l_upgo.ora.pdf',
            width = 55)

### GSEA GO -----------
sle.v3hvl.gsego <- sle_v3hvl |>
  filter(p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(OrgDb = 'org.Hs.eg.db',
        keyType = 'SYMBOL',
        ont = 'BP',pvalueCutoff = 1)

sle.v3hvl.gsego@result |>
  filter(NES > 0) |>
  plot_enrichment(metric = NES, padj_thres = .1)

sle.v3hvl.gsego@result |>
  filter(str_detect(Description, 'oxygen|stress|ndoplas')) |>
  as_tibble() |>
  select(1,2,NES,qvalue) |>
  DT::datatable()

sle.v3hvl.gsego |>
  gseaplot2('GO:0072593',
            title = 'GO pathway: reactive oxygen species metabolic process')

sle.v3hvl.gsego |>
  gseaplot2('GO:0034976',
            title = 'GO pathway: response to endoplasmic reticulum stress')

# v3h: SLE vs HC ---------
## FIG: volcano ===========
v3h_slevhc |>
  filter(p_val_adj > 0) |>
  plot_bill_volc(group1 = 'SLE V3-hi KC', group2 = 'HC V3-hi KC',
                highlights = c(key_cytokine, kegg_mva))

v3h_slevhc |>
  filter(p_val_adj > 0) |>
  plot_pub_volc(group1 = 'SLE V3-hi KC', group2 = 'HC V3-hi KC',
                highlights = c(key_cytokine, kegg_mva))

publish_pdf('mission/FPP/figures/v3h_SLE-HC_volcano.pdf')

## FIG: up ORA go ========
v3h_slevhc_upgo <- fine.type.slevhc.deg |>
  filter(p_val_adj < .05 & avg_log2FC > 1,
         str_detect(cluster, 'TRPV3-hi')) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Hs.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T,
           qvalueCutoff = .05) |>
  simpl_enrich()

v3h_slevhc_upgo@result |>
  publish_enrichment()

publish_pdf('mission/FPP/figures/mix4.epi.v3h.kc.SLE.vs.HC.upgo.ora.pdf',
            width = 60)

## GO GSEA -----------
v3h.slevhc.gsego <- v3h_slevhc |>
  filter(p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(OrgDb = 'org.Hs.eg.db',
        keyType = 'SYMBOL',
        ont = 'BP',pvalueCutoff = 1)

v3h.slevhc.gsego@result |>
  filter(str_detect(Description, 'oxygen|stress|ndoplas'), qvalue <.1) |>
  as_tibble() |>
  select(1,2,NES,qvalue) |>
  DT::datatable()

v3h.slevhc.gsego |>
  gseaplot2('GO:0072593',
            title = 'GO pathway: reactive oxygen species metabolic process\nV3h: SLE vs HC')

## in all skin cells ----------
epi_fine_fcs <-
  read_csv('mission/FPP/xiangya_sle_scRNA/epi.slevhc.all.logfc.csv')

### MVA logfc ---------
epi_fine_fcs |>
  filter(gene %in% kegg_mva) |>
  write_csv('mission/FPP/pub_source_data/SLE.vs.HC.skin.mva.logfc.csv')

epi_fine_fcs |>
  filter(gene %in% kegg_mva) |>
  ggplot(aes(gene, cluster, color = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point() +
  theme_pubr(x.text.angle = 45, legend = 'right') +
  scale_color_gradient2(low = 'blue', high = 'red')

epi.fine <- epi_fine_fcs$celltype |> unique()

skin.slevhc.gogse <- epi.fine |>
  map(\(x){epi_fine_fcs |>
      filter(celltype == x) |>
      pull(avg_log2FC, name = gene) |>
      sort(decreasing = T) |>
      gseGO(ont = 'ALL', OrgDb = 'org.Hs.eg.db',
            keyType = 'SYMBOL',pvalueCutoff = 1) |>
      pluck('result')}) |>
  set_names(epi.fine) |>
  list_rbind(names_to = 'cell.type') |>
  as_tibble()

skin.slevhc.gogse |>
  select(-c(setSize, enrichmentScore)) |>
  write_csv('mission/FPP/xiangya_sle_scRNA/results/skin.all.slevhc.go.gsea.csv')

skin.slevhc.gogse <-
  read_csv('mission/FPP/xiangya_sle_scRNA/results/skin.all.slevhc.go.gsea.csv')

skin.slevhc.gogse |>
  filter(str_detect(cell.type, 'hi'), qvalue < .05) |>
  write_csv('mission/FPP/xiangya_sle_scRNA/results/v3h.slevhc.gogse.csv')

plot.gsea.dot <- function(df) {
  df |> ggplot(aes(.data$cell.type, str_wrap(.data$Description, 50),
                   size = -log10(.data$pvalue), color = .data$NES)) +
    geom_point() +
    scale_size(range = c(0,3)) +
    scale_color_distiller(palette = 'RdYlBu') +
    theme_jpub +
    rotate_x_text(45) +
    labs(title = 'Differential pathway enrichment in SLE vs HC',
         y = 'GO term', x = 'Cell type')
}

skin.slevhc.gogse |>
  filter(str_detect(Description, 'reactive oxygen')) |>
  plot.gsea.dot()

publish_pdf('mission/FPP/figures/skin.slevhc.ros.gsea.pdf', width = 70)

skin.slevhc.gogse |>
  filter(str_detect(Description, 'TOR')) |>
  plot.gsea.dot()

publish_pdf('mission/FPP/figures/skin.slevhc.TOR.gsea.pdf', width = 70)

skin.slevhc.gogse |>
  filter(str_detect(Description, 'pyropto|inflammasome')) |>
  plot.gsea.dot() +
  scale_color_distiller(palette = 'RdYlBu', limits = c(-2.3,2.3))

publish_pdf('mission/FPP/figures/skin.slevhc.inflammasome.gsea.pdf', width = 75)

skin.slevhc.gogse |>
  filter(str_detect(Description, '[^o]autophag')) |>
  plot.gsea.dot()

publish_pdf('mission/FPP/figures/skin.slevhc.autophagy.gsea.pdf', width = 75)

skin.slevhc.gogse |>
  filter(str_detect(Description, 'unfolded protein')) |>
  plot.gsea.dot() +
  scale_color_distiller(palette = 'RdYlBu', limits = c(-1.7,1.75))

publish_pdf('mission/FPP/figures/skin.slevhc.UPR.gsea.pdf', width = 70)

skin.slevhc.gogse |>
  filter(str_detect(Description, 'toll')) |>
  plot.gsea.dot()

publish_pdf('mission/FPP/figures/skin.slevhc.TLR.gsea.pdf', width = 70)

skin.slevhc.gogse |>
  filter(str_detect(Description, 'interleukin-(1$|1 )')) |>
  plot.gsea.dot()

publish_pdf('mission/FPP/figures/skin.slevhc.IL1.gsea.pdf', width = 75)

## ERS-associated chaperone --------
erupr.list <-
  map_go_gene('GO:0034976',org = 'human')

unip.chap <-
  query_uniprot_keyword('KW-0143') |>
  separate_longer_delim(`Gene Names`, ' ') |>
  pull('Gene Names') |>
  unique()

erupr.chap <- erupr.list |>
  filter(SYMBOL %in% unip.chap) |>
  select(SYMBOL) |>
  write_csv('mission/FPP/UPR.associated.chaperone.gene.csv')

erupr.chap <-
  read_csv('mission/FPP/UPR.associated.chaperone.gene.csv')

fine.type.slevhc.deg |>
  filter(gene %in% erupr.chap$SYMBOL, str_detect(cluster, 'V3')) |>
  mutate(gene = fct_reorder(gene, p_val_adj, min),
         p_val_adj = ifelse(p_val_adj == 0, 1e-300, p_val_adj),
         cluster = fct_relevel(cluster, 'TRPV3-lo-KC')) |>
  ggplot(aes(y = cluster, x = gene, size = -log10(p_val_adj),
             fill = avg_log2FC)) +
  geom_point(shape = 21, stroke = .3) +
  scale_fill_gradient2(low = 'blue', high = 'red') +
  scale_size(range = c(0,3)) +
  theme_bw(base_size = 6) +
  labs(x = 'Gene', y = 'Cell type',
       title = 'ER stress-associated chaperon: SLE vs HC') +
  rotate_x_text(45)

publish_pdf('mission/FPP/figures/slevhc.V3KC.chaperone.logfc.pdf', width = 80,
            height = 40)

last_plot() |> pluck('data') |>
  arrange(gene, cluster) |>
  write_csv('mission/FPP/pub_source_data/sle.vs.hc.v3kc.chaperone.logfc.csv')

fine.type.slevhc.deg |>
  filter(gene %in% erupr.chap$SYMBOL, !str_detect(cluster, 'V3')) |>
  mutate(gene = fct_reorder(gene, p_val_adj, min),
         p_val_adj = ifelse(p_val_adj == 0, 1e-300, p_val_adj),
         cluster = fct_relevel(cluster, 'KC', after = Inf)) |>
  ggplot(aes(y = cluster, x = gene, size = -log10(p_val_adj),
             fill = avg_log2FC)) +
  geom_point(shape = 21, stroke = .3) +
  scale_fill_gradient2(low = 'blue', high = 'red') +
  scale_size(range = c(0,2)) +
  theme_bw(base_size = 6) +
  labs(x = 'Gene', y = 'Cell type',
       title = 'ER stress-associated chaperon: SLE vs HC') +
  rotate_x_text(45)

last_plot() |> pluck('data') |>
  arrange(gene, cluster) |>
  write_csv('mission/FPP/pub_source_data/sle.vs.hc.skin.chaperone.logfc.csv')

publish_pdf('mission/FPP/figures/slevhc.skin.chaperone.logfc.pdf', width = 80)

# V3l: SLE vs HC -----------
## FIG: volcano ==========
v3l_slevhc |>
  filter(p_val_adj > 0) |>
  plot_bill_volc(group1 = 'SLE V3-lo KC', group2 = 'HC V3-lo KC',
                 highlights = c(key_cytokine, kegg_mva))

v3l_slevhc |>
  filter(p_val_adj > 0) |>
  plot_pub_volc(group1 = 'SLE V3-lo KC', group2 = 'HC V3-lo KC',
                 highlights = c(key_cytokine, kegg_mva))

publish_pdf('mission/FPP/figures/v3l_SLEvHC_volcano.pdf')

## FIG: up ORA go ========
v3l_slevhc_upgo <- v3l_slevhc |>
  filter(p_val_adj < .05 & avg_log2FC > 1) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Hs.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T,
           qvalueCutoff = .05) |>
  simpl_enrich()

v3l_slevhc_upgo@result |>
  plot_enrichment() +
  ggtitle('Upregulated GO BP pathway in SLE V3-lo KC vs HC V3-lo KC')

v3l_slevhc_upgo@result |>
  publish_enrichment() +
  ggtitle('Upregulated GO BP pathway in V3-lo KC: SLE vs HC')

publish_pdf('mission/FPP/figures/v3l_SLE-HC_upgo.pdf', width = 65)

### GSEA GO -----------
v3l.slevhc.gsego <- v3l_slevhc |>
  filter(p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(OrgDb = 'org.Hs.eg.db',
        keyType = 'SYMBOL',
        ont = 'BP',pvalueCutoff = 1)

v3l.slevhc.gsego@result |>
  filter(str_detect(Description, 'oxygen|stress|ndoplas'), qvalue <.1) |>
  as_tibble() |>
  select(1,2,NES,qvalue) |>
  DT::datatable()

v3l.slevhc.gsego |>
  gseaplot2('GO:0072593',
            title = 'GO pathway: reactive oxygen species metabolic process\nV3l: SLE vs HC')

### DC: SLE vs HC up ORA go =======
dc_sle_upgo <- fine.type.slevhc.deg |>
  filter(avg_log2FC > 1 & p_val_adj < .05, cluster == 'DC') |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Hs.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           qvalueCutoff = .05,
           readable = T) |>
  simpl_enrich()

dc_sle_upgo@result |>
  publish_enrichment(scale_breaks = c(1e-23,3e-23,5e-23))

publish_pdf('mission/FPP/figures/mix4.epi.dc.SLE.vs.HC.upgo.ora.pdf',
            width = 55)

# LISA: UPR TFs ----------
upr_tf <- c('ATF4','ATF6','XBP1')

fine.type.slevhc.deg$cluster |>
  unique() |>
  walk(\(x)fine.type.slevhc.deg |>
        filter(p_val_adj < .05, avg_log2FC > 1, cluster == x) |>
        select(gene) |>
        write_source_csv(str_c('SLEvsHC_', x), col_names = F))

proj.nm <- 'mission/FPP/xiangya_sle_scRNA/'

lisa_slevhc <-
list.files('mission/FPP/xiangya_sle_scRNA/results/', 'lisa',
           full.names = T) |>
  read_tsv(id = 'file') |>
  mutate(file = str_remove_all(file, '.+HC_|.csv.+'))

lisa_slevhc |>
  filter(str_detect(factor, 'SREBF\\d$')) |>
  mutate(file = fct_relevel(file, 'KC', 'TRPV3-hi-KC', 'TRPV3-lo-KC')) |>
  ggplot(aes(file, -log10(summary_p_value), fill = factor)) +
  geom_col(position = 'dodge2') +
  labs(x = 'Cell type', fill = 'Transcription Factor',
       y = '-log10(p_value)',
       title = 'SREBF1/2 TF activity upregulated in SLE vs HC skin') +
  geom_hline(yintercept = -log10(.05), linetype = 'dashed') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub +
  rotate_x_text(45)

publish_source_plot('SREBFs.SLEvHC.skin.celltype.LISA', width = 80)

## HC all marker ---------
lisa_hc_all <-
  list.files('mission/FPP/xiangya_sle_scRNA/results/', 'motif_HC',
             full.names = T) |>
  read_tsv(id = 'file') |>
  mutate(file = str_remove_all(file, '.+marker_|.csv.+'))

lisa_hc_all |>
  filter(str_detect(factor, 'SREBF\\d$')) |>
  mutate(file = fct_relevel(file, 'Keratinocytes', 'TRPV3-hi-KC', 'TRPV3-lo-KC')) |>
  ggplot(aes(file, -log10(summary_p_value), fill = factor)) +
  geom_col(position = 'dodge2') +
  labs(x = 'Cell type', fill = 'Transcription Factor',
       y = '-log10(p_value)',
       title = 'UPR-associated TF activity upregulated in HC skin') +
  geom_hline(yintercept = -log10(.05), linetype = 'dashed') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub +
  rotate_x_text(45)

publish_source_plot('SREBFs.HC.skin.celltype.LISA', width = 80)

## SLE all marker ---------
lisa_sle_all <-
  list.files('mission/FPP/xiangya_sle_scRNA/results/', 'motif_SLE_all',
             full.names = T) |>
  read_tsv(id = 'file') |>
  mutate(file = str_remove_all(file, '.+marker_|.csv.+'))

lisa_sle_all |>
  filter(str_detect(factor, 'SREBF\\d$')) |>
  mutate(file = fct_relevel(file, 'Keratinocytes', 'TRPV3-hi-KC', 'TRPV3-lo-KC')) |>
  ggplot(aes(file, -log10(summary_p_value), fill = factor)) +
  geom_col(position = 'dodge2') +
  labs(x = 'Cell type', fill = 'Transcription Factor',
       y = '-log10(p_value)',
       title = 'UPR-associated TF activity upregulated in SLE skin') +
  geom_hline(yintercept = -log10(.05), linetype = 'dashed') +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  theme_jpub +
  rotate_x_text(45)

publish_source_plot('SREBFs.SLE.skin.celltype.LISA', width = 80)

